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Abstract We estimated a high-resolution local gravity field model over the south pole of the Moon 
using data from the Gravity Recovery and Interior Laboratory's extended mission. Our solution consists of 
adjustments with respect to a global model expressed in spherical harmonics. The adjustments are 
expressed as gridded gravity anomalies with a resolution of 1/6° by 1/6° (equivalent to that of a degree 
and order 1 080 model in spherical harmonics), covering a cap over the south pole with a radius of 40°. The 
gravity anomalies have been estimated from a short-arc analysis using only Ka-band range-rate (KBRR) 
data over the area of interest. We apply a neighbor-smoothing constraint to our solution. Our local 
model removes striping present in the global model; it reduces the misfit to the KBRR data and improves 
correlations with topography to higher degrees than current global models. 


Citation: 

Goossens, S., T. J. Sabaka, J. B. Nicholas, 
F. G. Lemoine, D. D. Rowlands, E. 
Mazarico, G. A. Neumann, D. E. Smith, 
and M. T. Zuber (2014), High-resolution 
local gravity model of the south pole 
of the Moon from GRAIL extended 
mission data, Geophys. Res. Lett., 41, 
doi:l 0.1 002/201 4GL0601 78. 


Received 9 APR 2014 
Accepted 8 MAY 2014 
Accepted article online 13 MAY 2014 


This is an open access article under 
the terms of the Creative Commons 
Attribution-NonCommercial-NoDerivs 
License, which permits use and distri- 
bution in any medium, provided the 
original work is properly cited, the 
use is non-commercial and no mod- 
ifications or adaptations are made. 


1. Introduction 

The Gravity Recovery and Interior Laboratory (GRAIL) mission was designed to map the structure of the lunar 
interior through high-precision global gravity mapping [ Zuber et al., 201 3a]. The mission was composed 
of two spacecraft with Ka-band intersatellite tracking as the single science instrument, complemented by 
tracking from Earth using the Deep Space Network (DSN) [Asmar etal., 2013]. The mission consisted of two 
phases: a primary mission which lasted from 1 March 201 2 to 29 May 201 2, during which the spacecraft had 
a mean altitude of 55 km above lunar surface, and an extended mission which lasted from 30 August 201 2 
to 14 December 2012. In the extended mission, the spacecraft had a mean altitude of 23 km above lunar sur- 
face until 1 8 November. After that, they were lowered to 20 km and finally to 1 1 km on 6 December. Initial 
analysis of the primary mission data resulted in a gravity field model in spherical harmonics of degree and 
order 420 [Zuber et al., 2013b], showing high correlations with topography up to high degrees and improv- 
ing on the pre-GRAIL lunar gravity models such as those by Konoplivet al. [2001], Matsumoto etal. [2010], 
and Goossens et al. [2011], which were based on Lunar Prospector, Selenological and Engineering Explorer 
(SELENE), and other historical data. Subsequent analysis of the primary mission data resulted in models of 
degree and order 660 [ Konopliv et al., 201 3; Lemoine et al., 2013], Owing to the much lower average alti- 
tudes above lunar surface in the extended mission, recent analysis has resulted in further improvements, in 
the form of models of degree and order 900 as presented by Konopliv etal. [2014] and Lemoine etal. [2014]. 
Both models have variable spatial resolution between degrees 600-680 over the central farside of the Moon 
to degree 900 at the poles. This variable spatial resolution of the gravity field is due to spacecraft altitude 
and ground track spacing [Konopliv et al., 2014]. Aliasing and relatively higher fits for the low-altitude part 
of the extended mission indicate that there is remaining signal in the Ka-band range-rate (KBRR) data even 
beyond degree 900. However, global spherical harmonics are not optimal when the spatial data coverage is 
heterogeneous, which is the case here at short wavelengths. 

In instances where global models require a large number of parameters to describe high-resolution features, 
the estimation of these parameters may become unstable, and local representations can be of great util- 
ity: they have the potential to extract high-resolution information from the data, while not demanding the 
computational burden of determining a full global model of similar resolution. In satellite geodesy there is a 
long history of local representations that have been applied successfully to different data types from satel- 
lites orbiting various planetary bodies. Examples include analysis for Venus [Barriotand Balmino, 1992], Mars 
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[Beuthe et al, 2006], and Jupiter's moon Ganymede [Anderson et a!., 2004; Palguta et a!., 2005], without even 
mentioning the many applications to Earth-orbiting satellites. Prior to the SELENE and GRAIL missions the 
tracking data distribution for the Moon was extremely asymmetrical with no farside data available, and thus, 
the Moon was especially suitable for the application of local methods, mostly using Lunar Prospector data: 
Sugano and Heki [2004a, 2004b] and Coossens et al [2005a, 2005b] estimated anomalies while Han [2008] 
and Han et al [201 1] used localized spherical harmonic functions to model the nearside gravity field. These 
same functions were used by Goossens et al. [201 2] to estimate the gravity field over the South Pole-Aitken 
basin on the farside of the Moon using SELENE data. In addition, local representations allow one to use con- 
straints other than the spectral constraints often used in global gravity modeling [e.g., Rowlands et al, 201 0; 
Sabaka et al, 2010]. 

The goal of the research presented here is to extract more information from the tracking data than can be 
modeled reasonably using global spherical harmonics. In the case of GRAIL KBRR data, our interest is espe- 
cially in areas where the altitude of the spacecraft was low and/or where the coverage was dense. Both 
Konopliv et al [2014] and Lemoine et al [2014] show that the degree strength of their models is at the max- 
imum degree at the poles because of the confluence of tracks. Therefore, we estimated a local gravity field 
model for the south pole of the Moon covering a cap of 40° radius, as adjustments to a global background 
model, at a resolution of 1/6° by 1/6° in both longitude and latitude, which is equivalent to the resolution of 
a degree and order 1 080 model in spherical harmonics (spatial block size = 5 km). 

2. Data and Methods 

GRAIL data consist of both DSN tracking (both in two-way and one-way modes) and intersatellite KBRR. 
While DSN data mostly determine the absolute positioning of the satellites, it is the contribution of KBRR 
data acquired at low altitudes that has resulted in the unprecedented GRAIL model sizes. KBRR data in the 
extended mission have a sampling of 2 s. Because we aim to estimate local parameters in one area of the 
Moon, we choose to use only KBRR data over the area of interest, excluding DSN tracking altogether. It is 
well known that each KBRR observation from the Gravity Recovery and Climate Experiment (GRACE) and 
GRAIL missions is approximately directly proportional to the difference in gravity potential at the two satel- 
lite locations. This makes KBRR data ideal for estimating local parameters from a subset of the entire mission 
data set [Rowlands et al, 2005]. KBRR data outside the region do not contribute, and data inside the region 
do not contain significant information from outside the region. 

Our analysis of the GRAIL extended mission KBRR data is based on the technique of Rowlands etal. [2002]. 
They showed that a short-arc analysis using KBRR data only is capable of recovering the underlying grav- 
ity field model. They used a baseline representation of the 1 2 parameters describing the difference vector 
between the two satellites and the state vector of the midpoint of the two satellites, instead of using the 
Cartesian inertial state vectors of the two satellites, all at the initial epoch. They showed that KBRR data are 
only sensitive to a subset of these baseline parameters, facilitating the coestimation of gravity parameters 
and orbit parameters from short arcs. This approach has been used successfully in the determination of 
monthly gravity solutions using GRACE data, for example Rowlands et al. [2005]; Luthcke et al. [2006], and 
later publications. 

Since we adjust gravity parameters on a spherical cap at the south pole, and since the GRAIL satellites were 
in polar orbits, each orbit passes over the area and each arc covering the area is approximately 25 min long. 
Using the orbits from a global gravity field model (GRGM900A, a predecessor to the GRGM900C model dis- 
cussed in Lemoine et al. [201 4]), we create short arcs over the south pole area with initial conditions and 
other parameters from these orbits and then adjust only the baseline parameters to the KBRR data. We 
emphasize that in this baseline analysis, we do not use the DSN data. The baseline parameters we adjust are 
those given in Rowlands etal. [2002] as being most sensitive to the KBRR data: the pitch of the GRAIL A-B 
position baseline, the magnitude of the GRAIL A-B velocity baseline, and the pitch of the GRAIL A-B veloc- 
ity baseline. We use the NASA Goddard Space Flight Center's (GSFC) GEODYN II Orbit Determination and 
Geodetic Parameter Estimation Program [Pavlis et al, 201 3] for the short-arc analysis. We adjust only the 
three baseline parameters for the short arcs, and once these arcs are converged we use GEODYN again to 
generate the partial derivatives of the KBRR data points with respect to the adjusted baseline parameters 
and the chosen gravity parameters. 
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We represent local gravity as gravity anomalies on a regular grid of 1/6° by 1/6°. This resolution is equivalent 
to that of a degree and order 1 080 expansion in spherical harmonics, and it was chosen in order to assess 
whether or not we could increase the resolution of our current global models (which are degree and order 
900). For a cap of 40° radius (with a 1/6° of a degree cap on the south pole itself), this leads to the estimation 
of a total of 518,401 anomalies over the selected area, for which we use 992,080 KBRR data points from the 
extended mission. This means that the ratio of observations and parameters is less than two and that con- 
straints will be important for the solution. We discuss the constraints below, and we discuss their influence 
on the solutions in section 3. Despite the fact that the longitudinal size of the grid cells physically decreases 
toward the pole, we did not compensate for this by decreasing the number of grid cells. Results in section 3 
will show we do not see artifacts in this dense south pole grid. 

Our gravity anomalies A g are related to the disturbing potential T (from W = U +T, where W is the gravity 
potential and U the normal potential of a reference field) in the following way: in a spherical approximation, 
the anomaly can be expressed as a function of T through A g = -dT /dr - 2 T /r, with r the spherical coordi- 
nate for radial distance. The potential T at a point P can be expressed in terms of anomalies A g(Q) at points 
Q on an exterior boundary <r as 


where R is the reference radius of the surface (1738 km in our case) and the kernel S(P, Q) is the Stokes-Pizetti 
kernel [ Heiskanen and Moritz, 1984]. If the disturbing potential is expressed in spherical harmonics, and the 
anomalies are global, then both anomalies and standard spherical harmonics are completely equivalent, 
and anomalies can be expressed in spherical harmonics, with the anomalies carrying a factor (/ - 1) in their 
expansion, where / is the degree. 

The partial derivatives of KBRR data with respect to the anomalies can be constructed from the accelera- 
tion on the satellite, by differentiation of equation (1). When computing the acceleration contribution per 
anomaly (from which the partial derivative is also computed), we approximate the integral in equation (1) 
by a summation. For the integration over each separate anomaly in equation (1), each block (of size 1/6° by 
1/6°) is divided into four subblocks (1/12° by 1/12°; except for 36 subblocks in the shape of wedges for the 
polar cap), and the integration is then approximated by the sum of the kernel function S(P, Q) evaluated 
at each center point of the subblock multiplied by the area of the subblock. This assumes that there are no 
anomaly contributions from outside of the chosen grid, which, in general, leads to boundary effects on the 
local solutions. This is further discussed in section 3. 

Once the partial derivatives of the data with respect to all parameters are generated, the normal equation 
system is formed. Due to the size of the problem, we use the supercomputers of the NASA Center for Climate 
Simulation (NCCS) at NASA GSFC. We apply a neighbor-smoothing constraint to our solutions following 
Rowlands et al. [2010] and Sabaka et al. [2010]. These spatial constraints are straightforward: each anomaly 
pair is forced to be equal, but each constraint equation is only given a weight W c) . that decreases with 
increasing distance between the /th and jth anomalies, according to W cjj = exp(1 - dJD ) with dy being 
the distance between the two anomalies (computed from their center points), and D being the correla- 
tion distance, which was here set equal to the latitudinal size of each grid cell (which is thus the same for 
all anomalies because the grid is regular). We discuss these constraints further in section 3. We apply the 
smoothing constraint to the full anomaly field, i.e., on Ag global + Ag adj where Ag g[obal is the vector with 
anomalies corresponding to the global spherical harmonics starting model and Ag adj is the vector with grav- 
ity anomalies for the adjustment. We apply the constraint to the full field (not just to the estimated gravity 
anomaly parameters) because the anomalies from GRGM900A show stripes (cf. Figure SI in the supporting 
information). Once the constraint system is formed, it is added to the data system with a scaling factor /<, 
resulting in the following equation: 


where A is the matrix with the partial derivatives of the data with respect to the parameters, W d is the data 
weight matrix, D is the matrix with partial derivatives of the constraints, W c the constraint weight matrix 
with elements W c y as defined above, r d is the vector with data residuals, and r c is the vector with con- 
straint residuals which is expressed as r c = -DAg g , oba |. The partial derivatives A and the data residuals r d 
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Figure 1 . Maps in stereographic projection centered at the south pole for (a and b) free-air anomalies (up to 50°) and 
(c and d) Bouguer disturbances (up to 70°S). The local model is shown in Figure la, the difference between the local 
model and GRGM900A is shown in Figure 1 b, GRGM900A (/ = 7-900) is shown in Figure 1 c, and a spherical harmonic 
transform of the local model (/ = 7-900) is shown in Figure Id. 

are evaluated using orbits determined with the global field GRGM900A. The anomalies for the global field 
are computed from its spherical harmonics expansion. The baseline parameters are not included in Ag adj 
but are accounted for by a modification of the data weight matrix W d (by taking into account their Schur 
complement), following Sabaka et al. [2010], The baseline parameters are unconstrained and free to adjust. 

3. Results 

We processed the GRAIL extended mission data as outlined in the previous section. We had 1 388 short arcs 
covering the selected south pole area. Our global starting gravity field model was GRGM900A. We weighted 
KBRR data at 0.1 pm/s throughout the entire extended mission because that is the level of fit with our cur- 
rent global models, but as in the processing of those models [Lemoine et al., 201 4], we noted improvements 
when the low-altitude data were relatively downweighted. A factor of 0.5 on the normal system for the 
December arcs was found to give the best results in terms of data fit, anomaly patterns, and correlations 
with topography (all of which will be shown and explained below). We also experimented with the weight 
of the constraint system, denoted as in equation (2), and we found that a value of n = 1 • 10~ 4 leads to the 
best overall results. We discuss this further at the end of this section. 

Figure 1 shows free-air anomalies and Bouguer disturbances (the latter of which are related to the dis- 
turbing potential through 8 g - —dT /dr, with a factor of (/ + 1 ) in their spherical harmonic expansion, as 
opposed tothe(/- 1) factor for anomalies) for our preferred local solution (Figures la and Id for free air and 
Bouguer, respectively), Bouguer disturbances for GRGM900A (Figure 1c), and anomaly differences between 
the local solution and GRGM900A (Figure 1 b). The free-air anomalies are shown up to 50°S, whereas the 
Bouguer disturbances are shown up to 70°S in order to focus on features close to the south pole. We also 
refer to Figures S2-S5 which show enlargements of these anomaly maps. The local model shows far fewer 
stripes than the global model, as the stripes in the difference map (Figure 1 b) negate stripes in the anomaly 
map of GRGM900A (cf. Figure SI ). We attribute this to our short-arc analysis where we adjust the orbit for 
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Figure 2. The root-mean-square (RMS) of fit to the KBRR data for GRGM900A and the local model, for the December arcs 
where the satellite altitude above topography was at its lowest. 

each track, and, more importantly, to the neighbor smoothing that was applied to the full anomaly field, 
thus including the anomalies from the starting model. We also tested solutions with either a damping con- 
straint (a diagonal constraint system, thus without neighbor smoothing) or with the neighbor smoothing 
applied to the adjustments only, and both solutions show more stripes in the anomalies. We also tested a 
neighbor-smoothing constraint with a varying correlation distance D (see section 2), where D was multiplied 
by the cosine of the latitude of the grid cell center point to account for varying size toward the pole. While 
this may slightly improve anomalies further away from the pole, the resulting small correlation distance at 
the pole leads to more striping in the solution there. An optimum between a constant and varying correla- 
tion distance may be found,but this was outside the scope of our current analysis. We also note that despite 
the decreasing physical size in longitude of the grid cells toward the pole, there are no apparent artifacts, 
likely thanks to the constant correlation distance applied. 

Bouguer disturbances were computed by subtracting the contributions of topography to gravity, expressed 
in the same principal axis frame as gravity, using the formulation of Wieczorek and Phillips [1 998]. We used 
Lunar Orbiter Laser Altimeter (LOLA) topography [Smith et at., 2010] archived in the Planetary Data System 
(LOLA/PDS Data Node, http://imbrium.mit.edu). We used a density of 2400 kg/m 3 to account for the 
observed decrease in crustal density with higher degrees [Wieczorek et at., 201 3; Han, 2013]. In order to com- 
pare two solutions of the same size so we can subtract the same topography-induced gravity potential, we 
generated global spherical harmonics for our local solution through a spherical harmonic transformation 
using Gauss-Legendre (GL) quadrature. We interpolated the adjustments Ag adj on a GL grid appropriate for a 
degree and order 1 079 model, applied the spherical harmonic transformation, and constructed an updated 
global model by adding the resulting coefficients to the background global model. Our local model thus 
becomes a global model with a maximum size of degree and order 1079 which is only attained in the south 
pole area; elsewhere it is equivalent to the background model. For the Bouguer disturbances in Figure 1 
we used both models up to degree 900. The Bouguer map for the local model shows fewer speckles and 
clearer features when compared to the map for the global model, which can help facilitate the geophysical 
interpretation of the structure of small-scale features [e.g., Zuber et at., 2014]. 

The free-air anomalies for the local solution in Figure 1 a show a nonphysical speckle pattern at longitude 
270°E, which is also where the adjustments are largest (Figure 1 b). This is the area just below Mare Orientale 
where GRAIL attained its lowest altitudes above topography. We have indications that for these areas the 
current resolution of 1/6° by 1/6° of a degree is not sufficient to extract the full signal in the data [Zuber et 
at., 201 4]. The speckle pattern is not as pronounced in the Bouguer disturbances in Figure 1, indicating that 
it appears beyond / = 900. 

The fits to the KBRR data for the global and local models, using short arcs over the selected area for both, are 
shown in Figure 2. The fits for the local model were computed by adding the acceleration contributions of 
the adjustments Ag adj using equation (1) to the global model GRGM900A in the short-arc orbit determina- 
tion. This shows that our local model fits the KBRR data better than the global starting model, sometimes by 
a considerable amount: on average, the fits improve by 1 2% for the whole extended mission and by 43% for 
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Figure 3. Localized correlations with topography-induced gravity for 
various gravity field models. The localization used a cap with a radius 
of 30° centered on the south pole. 


the low December arcs only. There are still 
high fits for certain periods, most notably 
around 5 and 8 December, where the 
GRAIL satellites had their lowest periap- 
sis heights. This indicates that information 
beyond degree and order 1080 is still 
present in the data, especially considering 
that the fit for those arcs is still well above 
the intrinsic noise levels of the data, which 
is thought to be between 0.05 and 0.07 
pm/s (N. Harvey, Jet Propulsion Laboratory, 
Pasadena, CA, unpublished data, 2013). 
Figure 2 shows the fits for the last part of 
the GRAIL extended mission because the 
improvements are most clear there. The 
fits also improve for the earlier arcs in the 
extended mission, although there they 
are very close to the global model fits, at 
— 0.12 pm/s. 


Figure 3 shows correlations of various models (we again used the spherical harmonic transform for the local 
solution) with the gravity potential generated from topography. We show correlations localized over the 
south pole, following Wieczorek and Simons [2005], using a cap of 30°, with one taper and a concentration 
factor of 0.9999, resulting in a windowing function with L win = 1 3. We chose the radius of 30° to avoid 
possible boundary effects at the edges of the solution, arising from the fact that we only use data over the 
area of interest, whereas data points close to the edge also contain information for unaccounted anomalies 
just beyond the grid. By testing small overlapping solutions, we found that the negative edge effects extend 
a few degrees at most, so our 1 0° cutoff is probably conservative. 


Figure 3 shows that the local model improves the correlations with topography considerably, and extends 
them to higher degrees. We also included the correlations for GRGM900C, a more recent model than 
GRGM900A, to show that the local solution also has higher correlations than that global model. Such high 
correlations are indicative of the pervasive fracturing of the lunar crust [Zuberet at., 2013b]. These results 
indicate that the smoothness of the maps in Figure 1 captures real gravity information and is not only the 
expression of the neighbor-smoothing constraint. 


Our solutions depend strongly on the data weights W d and the scaling factor p for the constraint system; 
cf. equation (2). We generated solutions with various weights for the December data and with various scal- 
ing factors //. We found that while anomaly maps for p = 1 • 1 0 -3 , in general, look smoother (they do not 
show as many stripes below the Mare Orientale area), such solutions do not fit the data well, nor do they 
maintain correlations with topography as high as the /i = 1 • 10~ 4 solution. To assess the contribution of 
the constraints to the solution, we decomposed the gradient component of equation (2) into two parts, one 
describing the data (A 7 W d r d ) and one describing the constraints (/<D r W c r c ). We then computed two sepa- 
rate solutions, each with only one of these components as a right-hand side to the equation system, which 
otherwise was the same for both solutions (thus including the constraint equation system /<D 7 W C D). These 
solutions indicate the contribution of each component to the final solution. These solutions are shown in 
the Figure S6, showing that the data largely influence anomalies in areas where the spacecraft were low 
and that the constraints largely influence the anomalies where GRGM900A showed stripes. Finally, down- 
weighting the December normal system reduces the striping in the area south of Mare Orientale, yet using 
a factor of 1 .0 rather than 0.5 slightly increases the correlations with topography. This again indicates infor- 
mation present in the data beyond our current resolution. However, the KBRR fits for such a solution in arcs 
other than the December arcs, especially earlier ones at slightly higher altitudes, are not as good as that for 
a solution using a factor 0.5, and often even higher than the fits for the starting model. For these reasons, 
a solution using a factor of 0.5 on the December system and a scaling factor of /< = 1 • 1 0^ 4 is our pre- 
ferred solution. For future solutions we might increase the grid resolution (in certain areas) to account for 
the information present in the low-altitude data. 
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4. Summary 

We presented a local model of the gravity field of the Moon over the south pole estimated from GRAIL 
extended mission data. Our model is expressed as gridded gravity anomalies with a resolution of 1/6° by 
1/6°, covering a cap over the south pole with a radius of 40°. We estimated the anomalies from a short-arc 
analysis using only KBRR data. We applied a neighbor-smoothing constraint to obtain our solution. Our local 
model removes striping present in the global models, and it improves fits to the KBRR data. We also trans- 
formed our local anomalies into spherical harmonics using Gauss-Legendre quadrature. This facilitates the 
comparison of Bouguer disturbances between the local and global models, and it allows us to compute 
correlations with topography. Our local model shows improved correlations with topography up to higher 
degrees than the global models. The spherical harmonic transformation from the original anomaly grid also 
has the benefit that it is easier to produce gravity disturbances and to compute gradients. Due to increased 
resolution and extended high correlations with topography, maps of the lunar gravity field such as from our 
local south pole model have the potential to improve the geophysical interpretation of small-scale features 
on the Moon. 
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